#ifndef _DISCRETE_FIELD_SYSTEM_H_
#define	_DISCRETE_FIELD_SYSTEM_H_

#include <map>
#include <numeric>
#include <memory>
#include <valarray>

#include <boost/math/constants/constants.hpp>
#include <boost/numeric/ublas/matrix.hpp>

#include "sync_queue.h"

#include "gl/glut.h"

#include <CGAL/squared_distance_2.h>
#include <CGAL/centroid.h>
#include <CGAL/Aff_transformation_2.h>

#include <mkl_blas.h>

#include "../common/polytope.h"
#include "../common/potential_energy.h"
#include "distance_field_grid_2.h"

/*
** Energy system computed with discrete field
*/

namespace Field_packing
{
	using namespace Packing_common;
	using namespace boost::numeric::ublas;

	class Discrete_field_system : public System_potential_energy
	{
		typedef Packing_polygon::Point_2 Point_2;
		typedef Packing_polygon::Vector	Vector_2;
		typedef CGAL::Line_2<CGAL::Exact_predicates_inexact_constructions_kernel>	Line_2;
		typedef CGAL::Aff_transformation_2<CGAL::Exact_predicates_inexact_constructions_kernel> Transformation_2;		

	public:
		Discrete_field_system(const std::vector<Packing_polygon*>& polygon_set, const Container_box& bndbox, Sync_queue< std::vector<Parameter> >& q) ;
		int optimize(double *trans, bool apply_constraints = false) ;
		int steepest_descent_optimize(double *trans);
		int directional_descent_optimize(double *trans);

		double system_energy_at(const double *const x);

		void energy_gradient_at(const double *const x, double *gradient) ;

		// objective function
		double gravity_energy_at(const double *const x) { return 0.0; }
		double molecular_energy_at(const double *const x) ;
		
		~Discrete_field_system()  {  }


		// debug
		static int newpoint_callback(const int evalRequestCode, const int n, const int m, const int nnzJ, 
							const int nnzH,	const double * const x,	const double * const lambda, 
							double * const obj, double * const c, double * const objGrad, double * const jac,
							double * const hessian,	double * const hessVector, void * userParams);

		void draw_sample_pnts() const;

		// ONLY FOR DEBUG
		void plot_graph();

	private:
		void constraint_at(const double *const x, double *constraint_values) ;
		void constraint_grad_at(const double *const x, double *constraint_grad) ;

		//void energy_gradient_at(const double *const x, double *gradient) ;
		void molecular_energy_gradient_at(const double *const x, double *gradient) ;

		void move_to_new_point(const double *const x);

	private:
		const std::vector<Packing_polygon*>& polygons;
		double rm; // at which the potential reaches minimum

		typedef std::vector<Point_2> Sample_set;
		std::vector<Sample_set> sample_pnts;

		// discretization grid
		std::auto_ptr<Distance_field_grid_2> dfg;
		//Distance_field_grid_2 *dfg;
		std::auto_ptr<Knitro_context> ktr_context;

		std::vector< matrix<double> > prev_trans;
		// parameters used in the attenuation function

		// debug
		//std::ofstream ofs;
		Sync_queue< std::vector<Parameter> >& trans_queue;
	};
}

#endif